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GEOMETRIC SOURCE SEPARATION SIGNAL PROCESSING TECHNIQUE 
1. FIELD OF THE IN VENTION 

This invention pertains generally to signal processing and more specifically to a 
system for performing mixed signal separation using geometric information and adaptive 
beamforming techniques. 



2. BACKGROUND OF THE INVENTION 

Blind source separation refers to the process of separating a composite signal into its 
15 original component signals without prior knowledge of their characteristics. This process is 
useful in speech recognition, multipath chamiel identification and equalization, improvement 
of the signal to interference ratio (SIR) of acoustic recordings, in surveillance applications and 
in the operation of hearing aids. 

20 Blind source separation of broad band signals in a multipatii environment remains a 

difficult problem which has a number of ambiguities. Increasing the number of sensors 
allows improved performance but leads to ambiguities in the choice of the separating filters. 
There are in theory multiple filters that invert the responses in a room because there are 
multiple projections fi-om the space containing microphone signals into the smaller space of 

25 signal sources. These multiple filters represent remaining degrees of fireedom in terms of a 

sensor array response. 



The consistent assigmnent of signal contributions to different source channels across 
30 different firequencies areates a fijequency permutation problem. This problem is inherent to j 
source separation algorithms including time domain algorithms unless the algorithm 
simultaneously considers different firequency bands. An estimation of such poiyspectral 
properties is particularly difficult for nonstationary signals such as speech and the resulting 
algorithms are computationally expensive. 
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The basic source separation problem is described by assuming the existence of M 
uncorrelated, time varying source signals source signals: 

s(t)€R^ 

5 

where the sources s(t) originate from different spatial locations. A number of sensors N 
(where N > M) detect time varying signals: 

x(t) € R^. 

10 

In a multipath environment each source j couples with sensor i through a linear transfer 
ftinction Ay(t), representing the impulse response of the corresponding source to sensor path 
such that: 

15 ^/fO=E>i Z^l'd Ay<'^)sx^-^)- 

This equation can be rewritten using matrix notation (denoting the convolutions by *): 
x(0-A(0*s(0. 

20 



After applying the discrete time Fourier transform, this equation may be rewritten as: 
25 = A(a))s(fi)). 

The goal of convolutive source separation is to find finite impulse response (FIR) filters 
W(/(t) that invert the effect of the convolutive mixing A(t). This is equivalent to producing 
y(fi?) = W(©)x(fl?) 

30 

that correspond to the original sources s(0- 
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Different criteria for convolutive separation have been proposed, for example as 
discussed by H.-L. N. Thi and C. Jutten in **BLIND SOURCE SEPARATION FOR 
CONVOLUTIVE MIXTURES", published in Signal Processing, vol. 45, no. 2, pp. 209-229 
(1995). A two channel example is disclosed in U.S. Patent No. 5,208,786 entitled MULTI- 

5 CHANNEL SIGNAL SEPARATION, issued to Weinstein et al. The '786 patent models each 
channel as a multi-input-multi-output (MIMO) linear time invariant system. The input source 
signals are separated and recovered by requiring that the reconstructed signals be statistically 
uncorrelated. However, the decorrelation condition is insxifficient to xmiquely solve the 
problem unless one assumes lhat the unknown channel is a 2 X 2 MIMO finite imptilse 

10 response filter. 

All convolutive separation criteria can be derived fi-om the assumption of statistical 
independence of the unknown signals, usually limited to pairwise independence of the source 
signals. Pairwise independence implies that all cross-moments can be factored, thereby 
15 creating a number of necessary conditions for the model signal sources. 



Ely^it)y']it^T)] = E[yUO]E\yJ(t^z)]. (EquationI) 

20 Convolutive separation requires that these conditions by satisfied for multiple delays x which 
correspond to the delays of the filter taps of W{x). For stationary signals higher order criteria 
(multiple n, m) are required. For non-stationary signals such as speech multiple ^ can be used 
and multiple decorrelation (n = m = 1) is sufficient. 

25 When using an independence criterion there remains both a permutation and scaling 

ambiguity. In the convolutive case the scaling ambiguity applies to each frequency group or 
bin resulting m a convolutive ambiguity for each source signal in the tkne domain. Any 
delayed or convolved versions of independent signals remain independent. For tiie 
independent frequency domain 

30 
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there is a permutation ambiguity per each frequency for all orders n and m. For each 
frequency the independent frequency domain is therefore also satisfied by arbitrary scaling 
and assignment of indices z, j to the model sources 

5 W{a))Ai(0)^P{6))Sia)), (Equation II) 

where P{co) represents an arbitrary permutation matrix and S(g)) an arbitrary diagonal scaling 
matrix for each frequency. This creates the problem that contributions of a given signal 
source may not be assigned consistently to a single model source for different frequency bins. 
10 A given model source will therefore have contributions from different actual sources. The 
problem is more severe with an increasing number of channels as the possible number of 
permutations increases. 

This problem has often been considered an artifact of the frequency domain 
15 formulation of the separation criteria since the separation task is decoupled into independent 
separation tasks per frequency bin. Forn = m = 1 this ambiguity also applies to the time 
domain independence criteria set forth in Equation 1. Even for higher orders the time domain 
criteria does not guarantee correct permutations. 

20 Some source separation work in the past simply ignored the problem. Others have 

proposed a number of solutions such as continuity in the spectra of the model sources, or the 
fact that the different frequency bins are often co-modulated. A rigorous way of capturing 
these statistical properties of multiple frequency contributions are polyspectra. 
However, in practice it is diflScult to obtain robust statistics at multiple frequencies, in 

25 particular for non-stationary signals such as speech. In addition, the algorithms that consider 
combinations of frequencies are inherently computationally very demanding. Smoothness 
constraints on the filter coefiBcients in the frequency domam have also been proposed, as for 

1 ^ T T o D«+*.r.* XTr» Alf^lAM PnfiflpH rONVOT J mVF. T^LTNm SOURCE 

SEPARATION USING A MULTIPLE DECORKELATION METHOD, issued to Parra et al. 

30 This is equivalent to constraining the length of the filter as compared to the size of the 

analysis window. However, this Umitation on the filter size may not always be reasonable as 
rather long filters are required in strongly reverberant environments. 
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In theory only N sensors are needed to separate M = N sources. In practice, however, 
one may want to use more microphones (N>M) to improve the performance of a real system. 
Ignoring the permutation and scaling ambiguities, Equation II reads 
W{a))A(co) = I, where /represents the identity matrix. For a given A^cd) there is a iV- M 
dimensional linear space of solutions Wicd), indicating that there are additional degrees of 
freedom when shaping a beam pattern represented by the jSlters W(ai). 

In conventional geometric and adaptive beamforming information such as microphone 
position and source location is often utilized. Geometric assumptions can be incorporated 
and implemented as linear constraints to the filter coefficients. In a multiple sidelobe 
canceler, for example, the response of one of the channels (channel i) is kept constant, which 
can be expressed as w(fi?)e ^ = constant The elements of the row vector w(fl?) e are the 
filter elements to be applied to each microphone, and 6/ is the iih column of the identity 
matrix. This is similar to the normalization condition imposed on the diagonal terms of W 
that is conventionally applied in blind separation algorithms. Rather than constraining a 
channel one can also constrain the response of a beamformer for a particular orientation. 

If the locations and response characteristics of each microphone is known, one can 
compute the free field response of a set of microphones and associated beam forming filters 
w(fl>). For a position q, the phase and magnitude response is given by 

r(a>,q) = w(fi>)d(<»,q), 

where d(a?,q) e represents the phase and magnitude response of the iV microphones for a 
source located at q. For a linear array with omnidirectional microphones and a fer field 
source (much beyond the array aperture squared over the wavelength of interest) the 
microphone response depends approximately only on the angle G = e(q) between the source 
and the linear array. 
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where pi is the position of the ith microphone on the linear array and c is the wave 
propagation speed. 

Constrauiing the response to a particular orientation is sinq)ly expressed by the linear 
constraint on w(<») such that r(fi>,9) = w(fi))d(fl?,e) = constant. This concept is used in the 
linearly constiained minimum variance (LCMV) algorithm and is also the underlying idea for 
generalized sidelobe cancelling. In order to obtain a robust beam it has also been suggested to 
require a smooth response around a desired orientation. In summary, all of these conditions 
or a combination of tiiem can be expressed as linear constraints on w(o). 

Most adaptive beamforming algorithms consider power as tiieir main criteria for 
optimization. Sometimes power is minimized such as in noise or sidelobe cancelling in order 
to adaptively mimmize tiie response at the orientation of tixe interfering signals. Sometimes 
power is maximized, as in matched filtw approaches, to maximize the response of interest 
As a result, tiiese algorithms often perform suboptimally when there is cross talk from other 
sources. 

In second order soiurce separation methods rather than considering the power of an 
individual beam w(fi)) e C "'^and an individual channel y(t) e R', one can consider powers 
and cross powers of multiple beams W(a)) e C '^''^and tiieir corresponding outputs y(t) e R 
. In flie frequency domain these multiple beams and outputs correspond to the cross power 
spectra R3^T,a>) . Second order blind source separation of nonstationary signals minimizes 
cross powers across multiple times. Off diagonal elements of tiie matrix R»<T,fi>j are 
minimized in second order separation rather tiian flie diagonal terms as is the case in 
conventional adaptive beamforming. Stict one chamiel power criteria has a serious cross talk 
or leakage problem when multiple sources are simultaneously active, especially in reverberant 
environments. 
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SUMMARY OF THE INVENTION 

The present invention addresses the need to separate the source of nonstationary broad 
band signals in a multipath environment. Ambiguities in source separation are addressed by 
adding prior information such as microphone position and by adding the assumption that the 
sources are localized in space. By replacing the power criteria of many adaptive 
beamforming algorithms by a cross power criteria a plurality of geometric source separation 
algorithms are obtained. To address the permutation, convolution and more general degrees 
offreedom forthe shaping ofbeams, geometric information is utilized. The advantages of 
bllind source separation and geometric beamforming are combined by minimizing the cross 
power spectra for multiple times t while enforcing constraints used in conventional adaptive 
beamforming. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Figure 1 depicts a system for executing a software implementation of the present 
invention; 

Figure 2 is a flow diagram of a method of the present invention; 

DETAILED DESCRIPTION OF THE INVENTION 

The present invention estimates values of W{ai) obtained by known blind source 
separation techniques by making geometric assumptions regarding the signal sources. The 
sources are assiuned to be localized at least up to the spatial resolution of a given array. The 
invention assumes that signals originate from the same locations for the entire frequency 
spectrum, permitting the formulation of geometric constraints on the filter coefficients. 
Dififerent geometric constraints are introduced leading to a class of Geometric Source 
Separation algorithms. 

Figure 1 depicts a system 100 for implementing the source separation method of the 
present invention. The system 100 comprises a composite signal source 126 that supplies the 
signal that is to be separated into its component signals and a computer system 108 that 
executes the geometric source separation routine of the present invention. The source 126 
may contain any source of composite signals, but is illustratively shown to contain a sensor 
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array 102, a signal processor 104 and a recorded signal source 106. The sensor array contains 
one or more transducers 102A, 102B and 102C such as microphones. The transducers are 
coupled to a signal processor 104 that performs signal digitization. A digital signal is coupled 
to the computer system 108 for signal separation and fiulher processing. A recorded signal 
5 source 106 may optionally form a source of the compositie signals that require separation. 

The computer system 108 comprises a central processing unit (CPU) 1 14, a memory 
122, support circuits 1 16 and an input/output (I/O) interface 120. The computer system 108 is 
generally coupled through the I/O interface 120 to a display 1 12 and various input devices 
10 110 such as a mouse and keyboard. The support circuits generally contain well known 

circuits such as cache, power supplies, clock circuits, a communications bus and the like. The 
memory 122 may include random access memory (RAM), read only memory (ROM), disk 
drive, tape drive and the like, or some combination of memory devices. 

15 The invention is implemented as the geometric source separation routine 124 that is 

stored in memory 122 and executed by the CPU 1 14 to process the signal from the signal 
sources 126. As such the computer system 108 is a general purpose computer system tiiat 
becomes a specific purpose computer system when executing the.routine 124 of the present 
mvention. Although a general purpose computer is illustratively shown as a platform for 

20 implementing the invention, the invention can also be implemented in hardware as an 
application specific integrated circuit (ASIC), a digital signal processing (DSP) integrated 
circuit, or other hardware device or devices. As such, the invention may be implemented in 
software, hardware or a combination of software and hardware. 

25 Figure 2 depicts a flow diagram of the geometric source separation routine 124 of the 

present invention. At step 200, the composite, mixed signal 126 is input, the signal being 
separated into a plurality of data frames contaming data san5)les of the input signal x(t). The 
routine 200 produces discrete Fourier transform (DFT) values x(ffl) for each data frame jc(t), 
that is, one frequency domain DFT value for each wmdow of length T san^les. 

30 

At step 204, a running time estimate of Ryy(«,ffl) is computed from the outputs y(0- As 
large T are required for typical filter sizes, simultaneous diagonalization of Ryy(r,©) for 
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multiple times t is perfbmied so as to minimize the simi of squares of the off diagonal 
elements according to the equation 213: 

AW)= 2 Ryy(^^^)-diag(Ryyaa>)) 1 ' , 

5 

with the factorization Ryy(/,6)) « W(fi)) Rxx(^£y) W"(fl>). || ... || refers to the Frobenius norm 
defined as || M || ^ = Tr (MM^). The siramiations over t and <d will range over all time and 
all frequency bins for which adaptation of W will occur, respectively. For faster 
convergence using gradient descent the total input power per frequency a {oS) is normalized, 

10 

This criteria is minimized with respect to the filter coefficients W. The lower bound of zero is 
obtained only if Ryy(^ty) is diagonal. 

15 

The signal sources s(0 are localized at angles 0 = [Gi ... Gm] and at a sufficient distance 
from the sensors 102 for a far field approximation to apply. While full three dimensional 
source location can be advantageously used with the present invention, the specific 
embodiment described here will identify souce locations simply by incident angle to the 
20 microphone array. Step 208 computes the geometric source locations based on the cross 
power spectra produced at step 204. Step 219 applies a selection or switch to determine 
which of various criteria will be next be applied to the cross power spectra calculation. For 
example, if the hard constraint 210 



25 diag iWiai) D(a), G)) = 1 

is applied to each filter W/(fi)) , then flie rth row vector in W(a)) is forced to have a unit 
response in tiie direction G,. In a gradient descent algorithm constraint 210 can simply be 
implemented by projecting the gradient {SJIdH^iico) ) to the Unear subspace of permissible 
30 solutions with a constrained gradient. Smce power or cross power mimmization will try to 
minimize the response at the interference angles, this will lead to an equivalent singularity at 
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20 



30 



10 

those frequencies. In those cases soft constraint 210 should be selected as a regularization 
term 211 of the form 

/,(fi,)=| diag(in®)2)(t»,e))-l 

A second, more resrictive hard the constraint 209 
FP(fi))D(o,e) = / 

10 may be selected a step 219. Hard constraint 209 imposes the conditions of hard constraint 
210 to each filter w<o) by forcing the ith row vector in Wito) to have a unit response in the 
direction 6,. Further, the hard constraint 209 requires that the ith row vector have a zero 
response in the direction of interfering signals 6/ where i ^j. 

15 At those frequencies where the grating lobes (periodic replica of the main lobes due to 

limited spatial sampling) of a beam pattern cross the interfering angles, D(<», 6) is not 
invertible. In those cases it is unreasonable to try to enforce constraint 209 as a hard 
constraint. Rather, soft constraint 212 is selected by adding the regularization term of flie 
form 



Step 219 performs various intialization conditions on equation 213. In all initializations the 
relationship Wi(<o)e/ = 1 is required to normalize the scale during optiinization. In step 215the 
25 filter structure is initialized to correspond to a delay-sum beamformer pointing into the 
orientations of the individual sources (del-sum). Using the orientations Of the filter 
«««.«;^;<.„tr «f tv.A ith h^am fthf. mw vRTitfirs W of W( mW are initialized with algorithm 



zeroes 



Wi(G>)= d(a),e/ . 

In step 219 flie equation 213 can alternatively be initialized witii beams that place 
at all orientations of interfering sources, that is, at angles B/Gifor tiie ith beam. The 
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initialization filters with minimum nomi tbat satisify those conditions can be computed 
explicitly with a least squares approach resulting in the mitialization algorithm 216, 

5 . ' • ' ■ , " • . 

where indicates the psuedo-inverse and [e,-, D{a), 0/0,)] is a matrix containing the zth column 

6/ of the unit matrix and all but the ith column of D(fi?, 0). 

In an on line implementation of a separation algorithm the concept of introducing 
10 geometric information through an initialization is tj^ically not feasible as the source positions 
in the environment may be changing dynamically. More frequently the geometric constraints 
are enforced at initialization and throughout the optimization process. The constraints change 
dynamically as the estimated locations change. The linear constraints 209 and 210 are each 
typically implemented as a soft constraint with a penalty term. The problem of 
15 noninvertibility is addressed by introducing a frequency depeiident weighting of the penalty 
term. In particular the goal is to eliminate the constraints from the optimization for those 
frequency bands for which 2>(fi>, 0) is not inveritble. A rather straightforward metric for 
invertibiUty is the condition number. Therefore at steps 211 and 212 the regularization term 
J(W) is weighted with the the inverse of the condition number of T^iai) = cond" 
20 0)), which converges to zero when D(a), 0) is not invertible and remains bounded 

otherwise, for example as 0 < X(a)) < 1 . The total cost function 218 including frequency 
dependent weighting of the geometric regularization term is given by 

25. 

In algorithm 210 the regularization term J\ will try to maintain the response of filters / in 
orientation 0,. The delay-sum beamfonner 215 satisfies the condition of algorithm 210 
strictly. In algorithm 209 the regularization term J2 will ia addition minimize the response for 
the orientations of the interfering sources. The filter stracture that strictiy guarantees the 
30 constraints of algorithm 209 are computed with a least squares approach as the psuedo inverse 
of D^(g>, 0), or by including a regularization term /ST for the noninvertibility problem, 
resulting in algorithm 219: 
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W{(0) ^D\co, 9) 9) D\(D, 9) + pi)'^' 

Algorithm 2 1 7 places a zero at the angles of interfering sources, while its response in 
5 other directions is not specified. The results for algorithms 2 1 5, 210 and 209 exhibit a main 
lobe in the directions of the corresponding source. For conflicting frecjuency bands, where a 
grating lobe coincides with the location of an mterfering source, multiple cross power 
minimization cancels the main lobe for algorithm 215, while conserving it somewhat for 
algorithms 210 land 209 due to the geometric penalty. Qualitatively, the results for the data 
10 independent algorithm 219 appears to capture both main lobe and zeros in the correct 

location. However, its performance is inferior to the data adaptive algorithms 209, 210, 215 
and 217. 

Any of the algorithms 209, 210, 21 1 and 212 can be used to optimize the signal filter 
15 222 and minimize cross power at step 202 using a gradient descent algorithm.. The angles 9, 
of the multiple sources can be identified automatically using the Multiple Signal 
Classification (MUSIC) algorithm, which is a metiiod of computing directions of arrival of 
multiple signals using arbitrarily placed antennas with arbitrary directional responses. The 
MUSIC algorithm also applies to measuring firequencies of multiple sinusoids comprising a 
20 sampled data time series. When used to measure frequency it is directly suited to tapped 
delay line implementations. The number of sources M is assumed to be known. While the 
examples and algorithms given in this description are directed toward numerical operations in 
the frequency domain, the preset invention can also be implemented in the time domain. 

25 



30 
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CLAIMS 

L An apparatus for separating a convolutively mixed signal into a plurality of sources, 
comprising: 

5 a plurality of signal sensors ( 1 02), each sensor residing at a sensor location; 

a multiple input, multiple output signal filter (222); 

means (204) for estimating cross powers of multiple output chamiels at multiple times; 
means (208) for computing a desired set of source locations r, the source locations r 
being defined by spatial relation to the sensor locations; and 
10 means for adapting the signal filter to minimize midtiple output powers while 

enforcing a desired filter response for the desired set of source locations. 

2. The apparatus of claim 1 wherein the cross powers are minimized by estimating and 
minimizing multiple cross power spectra of the signal filter outputs in a fi'equency domain. 

3. The apparatus of claim 1 wherein the desired filter response is enforced as a linear 
constraint on the signal filter when the signal filter is minimizing multiple output powers. 

4. The apparatus of claim 1 wherein the desired filter response is enforced by adding a 
20 regularization tenh when the signal filter is minimizing multiple output powers, 

5. The apparatus of claim 2 wherein a imit response for source locations r are enforced as 

diag(FF(fi>)2)(fi>,r))=l 

25 6. The apparatus of claim 2 wherein a xmit response and a zero response for soiurce locations r 
are enforced as 

FF(fi>)D(G3!,r) = I 

7. The apparatus of claim 1 wherein the signal filter (222) is a finite impulse response filter. 
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8. The apparatus of claim 2 wherein cross power minimization is implemented by minimizing 
a cost function (213) defined as 



5 



J(W)= £ ^(^) II Ryy(^,^)"<iiag(Ryy(^^i>))|| 



9. The apparatus of claim 2 wherein a unit response for source locations r is enforced with an 
additional regularization term (211) defined as 



10. The apparatus of claim 2 wherein a imit response and a zero response for locations r are 
enforced with an additional regularization term (212) defined as 



1 1 . The apparatus of claim 1 wherein the source locations r are specified by angles relative 
to a signal sensor location. 

20 12. The apparatus of claim 1 wherein the signal filter (222) is initialized to behave as a beam 
former filter. 

13. The apparatus of claim 12 wherein the initialized signal filter represents a set of delay 
sum beamformers aiming at the desired set of source locations r. 



14. The apparatus of claim 1 wherein the source locations r are estimated firom cross power 
spectra. 



10 



/i(W)= II diag(W(fi))2)(fl>,r))-l|| ' . 



15 



/2(W)=|| (r(fi?)2>(fi?,r)).l|| ' . 



25 
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15 . A computer system (1 08) for processing nonstationaiy signals, comprising: 

at least one signal input port; 

a central processing unit (1 14) configured to receive a composite signal from the 
signal input port; 

a memory (122) interconnected to the central processing unit (1 14); and 
a geometric source separation modvUe (124) interconnected to the memory (122) to 
process digitized signals stored in the memory. 

16. The computer system of claim 15, wherein the geometric source separation module 

jRirther comprises at least one finite impulse response filter (222), the composite 
signal being applied to the filter. 

17. The computer system of claim 16, wherein the memory (122) further comprises 

spatial location data of sources of the composite signal, said spatial location data 
being coupled to the finite impulse response fiilter (222). 

18. The computer system of claim 17, wherein the geometric source separation module 

furflier comprises a pluraUty of algorithms (209, 210, 215, 217) adapted to 
compute a plurality of filter coefficients for the finite impulse response filter. 

19. A method of separating a composite signal into a plurality of component signals 

comprising the steps of: 

computing the cross power spectra Kyyit,co) of the plurality of component signals; 
simultaneously diagonalizdng Ryy(<,ffl>) with an algorithm (213); 
applying a linear constraint (209,210) to the algorithm (213) to create filter 
cuefQcients; and 

filtering the composite signal based on the filter coefficients. 
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16 



20. The method of claim 19, flirther comprising the step of defining the diagonalizing 
algorithm (213) as 

J(W)= £ a(i») I Ryy(f,©)-diag(Ryy(/,a?)) II ' , 
5 with the factorization RyyM » .W(fi?) R^xM W"(fi?). 

21. The method of claim 20, further comprising the step of choosing a linear constraint 

from a set of two linear constramts given by a first linear constraint (210) 
diag(»f{(»)D(fi),e)) = l 
10 and a second linear constraint (209) 



22. The method of claim 21, further comprising the step of applying a regualarization term 
(21 1) to the diagonalizing algorithm (213) of 
15 m = II diag (W(0)D{G>, 9)) - 1 P . 
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